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ABSTRACT 

1) Context. Protostellar accretion discs have cool, dense midplanes where externally originating ionisation sources such as X-rays 
or cosmic rays are unable to penetrate. This suggests that for a wide range of radii, MHD turbulence can only be sustained in 
the surface layers where the ionisation fraction is sufficiently high. A dead zone is expected to exist near the midplane, such 
£Nj that active accretion only occurs near the upper and lower disc surfaces. Recent work, however, suggests that under suitable 
conditions the dead zone may be enlivened by turbulent transport of ions from the surface layers into the dense interior. 
Aims. In this paper we present a suite of simulations that examine where, and under which conditions, a dead zone can be 
enlivened by turbulent mixing. 

n Methods. We use three-dimensional, multifiuid shearing box MHD simulations, which include vertical stratification, ionisation 
O chemistry, ohmic resistivity, and ionisation due to X-rays from the central protostar. We compare the results of the MHD 
+3 simulations with a simple reaction-diffusion model. 

& Results. The simulations show that in the absence of gas-phase heavy metals, such as magnesium, turbulent mixing has essentially 

1 i no effect on the dead zone. The addition of a relatively low abundance of magnesium, however, increases the recombination 

time and allows turbulent mixing of ions to enliven the dead zone completely beyond a distance of 5 AU from the central star, 
t-H for our particular disc model. 
k> Conclusions. During the late stages of protoplanetary disc evolution, when small grains have been depleted and the disc surface 
density has decreased below its high initial value, the structure of the dead zone may be significantly altered by the action of 
turbulent transport. This may have important consequences for ongoing planet formation in these discs. 
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1. Introduction may prevent settling toward the midplane (Carballido et 

al. 2005; Johansen & Klahr 2005; Fromang & Papaloizou 
Observations of young stars in a broad variety of star 2 006). It will also cause an increase in planetesimal veloc- 
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forming environments show that they are surrounded by ity dispersion, and may prevent runaway growth of plan- 

Vh gaseous and dusty circumstellar discs (e.g. Beckwith & etesimals into planetary embryos (Nelson 2005). It will 

Sargeant 1996; O'Dcll et al. 1993; Prosser et al 1994; mo dify the migration of low mass protoplanets (Nelson & 

Furlan et al. 2006). These systems usually show signatures p apa i izou 2004; Nelson 2005), and will provide the effec- 

of active gas accretion onto the central star with a range of tive v j scous stre ss in the disc needed to drive the so-called 

mass flow rates, but with the typical value being ~ 10~ 8 typc n migration of giant planets (Nelson & Papaloizou 

M yr" 1 (e.g. Sicilia-Aguilar et al. 2004). The mechanism 2 003). It has recently been suggested that it may lead to 

by which the disc transports angular momentum inter- planetesimal formation through gravitational instability 

nally to cause accretion is not yet fully understood, but is (Johansen et al 2007) 
likely to have its origin in disc turbulence. So far only one 

mechanism has been shown to be robust in generating There are continuing questions, however, about the ap- 

turbulcnce in Keplerian discs, namely the magnetorota- plicability of the MRI to cool, dense protostellar discs, 

tion instability (MRI) (Balbus & Hawley 1991; Hawley & as the ionisation fraction is expected to be low (Blaes 

Balbus 1991). In addition to providing the internal stress & Balbus 1994). A protostellar disc model has been pre- 

rcquired for accretion, turbulence may also have impor- sented by Gammie (1996) in which the main source of 

tant consequences for planet formation and evolution in ionisation is Galactic cosmic rays. Such a disc is predicted 

protoplanetary discs. It will lead to mixing of the dust, and to have magnetically "active zones" near the disc surface 
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in which turbulence is sustained by the MRI due to cos- 
mic ray ionisation, but with a "dead zone" near the disc 
midplane where cosmic rays are unable to penetrate. Sano 
et al. (2000) have examined the effects of a more complex 
chemical reaction network and the influence of small dust 
grains. Glassgold et al. (1997) and Igea & Glassgold (1999) 
examined X-rays emitted by the protostellar corona as a 
possible source of disc ionisation, since it is doubtful that 
cosmic rays can penetrate the inner disc regions because 
of the attenuating effect of the T Tauri wind. Fromang et 
al. (2002) demonstrated the potential importance of gas 
phase heavy metals such as magnesium, whose presence 
in trace quantities can significantly increase the recom- 
bination time and decrease the size of the dead zone (at 
least in a dust free disc). Semenov et al. (2004) examined 
the disc chemistry and ionisation fraction using a reaction 
network drawn from the UMIST database, and analysed 
results using a reduced reaction network. 

Nonlinear numerical studies of the effects of ohmic re- 
sistivity on MRI-driven MHD turbulence have also been 
presented. Fleming et al. (2000) performed MHD simu- 
lations of resistive discs, and showed that turbulence is 
not sustained in discs with a magnetic Reynolds num- 
ber Re m which is below a critical value Re™* . Fleming & 
Stone (2003) performed simulations where resistivity de- 
creased as a function of height in the disc, as predicted by 
the Gammie (1996) model, and showed that active zones 
could indeed coexist with dead zones in the disc. They 
also showed that a low Reynolds stress can be maintained 
in the dead zone, such that low levels of accretion are sus- 
tained there. A recent study of resistive discs has also been 
presented by Turner et al. (2007), who presented stratified 
shearing box simulations of discs in which resistivity var- 
ied with height, and a multifluid simulation in which disc 
chemistry was evolved along with the dynamics. This lat- 
ter run showed that turbulent mixing can potentially have 
an important effect in generating stresses in the dead zone. 
Fromang & Papaloizou (2007) have recently performed a 
resolution study of shearing box simulations, and showed 
that the level of turbulent activity reduces as the resolu- 
tion increases. An analysis of exiting shearing box simu- 
lations in the literature by Pessah et al. (2007) led to a 
similar conclusion. In a companion paper to Fromang & 
Papaloizou (2007), Fromang et al. (2007) also examined 
how turbulent activity scales with magnetic Prandtl num- 
ber (defined by P m — Re m /Re where Re is the Reynolds 
number and Re m is the magnetic Reynolds number) . They 
showed that both the Reynolds number and magnetic 
Prandtl number are the parameters that control the level 
of turbulent activity in a disc, with P m < 1 flows showing 
no sustained turbulent activity for zero net flux magnetic 
fields. Clearly there is much work to be done in under- 
standing the nature of MHD turbulence in discs. 

In a recent set of publications, we have undertaken 
an extensive study of the chemistry and ionisation frac- 
tion in protoplanetary discs. In Ilgner & Nelson (2006a) 
(hereafter paper I) we examined the dead zone structure 
in standard a-disc models as predicted by a number of 



different chemical reaction networks, and showed that the 
simple model of Oppenheimer & Dalgarno (1974) gives 
good agreement with more complex models based on the 
UMIST database (Le Teuff et al. 2000). We also demon- 
strated that grain depletion by factors even lower than 
10~ 6 are required to reduce significantly the sizes of dead 
zones. In Ilgner & Nelson (2006b) (hereafter paper II) 
we constructed a reaction-diffusion model to examine the 
role of turbulent mixing on dead zone structure. The re- 
sults showed that turbulent mixing has essentially no ef- 
fect throughout the disc in the absence of gas phase heavy 
metals such as magnesium. In the presence of trace quan- 
tities of magnesium, however, turbulent mixing was able 
to enliven the dead zone out beyond a few AU, where the 
mixing time becomes shorter than the recombination time. 
In a third paper of the series (Ilgner & Nelson 2006c), we 
examined the effect of X-ray flares on dead zone structure. 

In this paper we present a suite of shearing box sim- 
ulations of stratified local disc models in which we evolve 
the disc chemistry along with the magnetohydrodynamic 
equations. The primary aim is to re-examine the results of 
paper II using multifluid MHD simulations, and determine 
whether, and under what conditions, turbulent mixing can 
enliven a dead zone by mixing ions from the surface layers 
down toward the midplane. We use the simple reaction 
scheme of Oppenheimer & Dalgarno (1974), which we in- 
corporate into a multifluid MHD code, and assume that 
dust grains are absent and ionisation is caused primarily 
by X-rays from the central star. We examine the effects 
of mixing as a function of gas phase magnesium abun- 
dance and distance from the central star. Our results are 
in very good agreement with the predictions of paper II. 
Disc models which contain no gas phase magnesium show 
that the dead zone structure is essentially unmodified by 
turbulent mixing. In the presence of magnesium, however, 
our simulations show that the dead zone can be enlivened 
completely beyond a distance of 5 AU from the central 
star. 

The paper is organised as follows. In Sect. [2] we present 
the basic equations and the chemical reaction network 
that we solve. In Sect. [3] we discuss the reaction-diffusion 
model, which we use to compare with the MHD simula- 
tions. In Sect.|4]we discuss the method used for calculating 
the X-ray ionisation rate, and in Sect. [5] we discuss previ- 
ous simulations that have examined dead zone structure. 
In Sect.[6]we present our simulation results, and finally in 
Sect. [7] we draw our conclusions. 

2. The dynamical and chemical model 

In this section we give a detailed description of the chem- 
ical model used in our simulations, and present the mul- 
tifluid MHD equations that we solve. 

2.1. Chemical model 

For the purposes of simplicity and computational 
tractability, we have applied the simple kinetic model 
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Table 1. Rate coefficients for the Oppenheimer & 
Dalgarno model. 
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of Oppenheimer & Dalgarno (1974) to evolve the gas- 
phase chemistry within the simulations. This reaction 
network has been described in Ilgner & Nelson (2006a), 
where it was compared with more complex reaction net- 
works and found to predict electron fractions that were 
slightly higher on average, due to the lower number of 
molecular ion species present in the simpler model. The 
Oppenheimer & Dalgarno reaction network may be writ- 
ten: 



H 2 + 
H+4 
H+4 
Mg+ 



hv 
e~ 

■Me 



H. 
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H 2 

H 2 - 
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Mg+ 
+ hv 
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(3) 
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This reaction scheme involves five species: molecular hy- 
drogen and its ionised counterpart (which act as a repre- 
sentative molecule and molecular ion) , atomic magnesium 
and its singly ionised counterpart (which act as represen- 
tative heavy metal and metal ion), and free electrons. The 
free electrons are treated as a dependent species and their 
local fractional abundance is calculated assuming local 
charge neutrality: x[e~] = x[H^] + a;[Mg + ]. In our model, 
Eqns. ([l]) — (j4j) represent ionisation of H 2 by X-rays, re- 
combination of , charge transfer between H^ and Mg, 
and recombination of Mg + , respectively. The associated 
reaction rates are given in table [I] with the rate coeffi- 
cients being denoted by £, a, (3, 7, respectively. The re- 
combination rate of Mg + is five orders of magnitude lower 
than of H 2 + , indicating that a magnesium-abundant gas 
will sustain a higher ionisation fraction than a metal-poor 
one because of charge transfer reactions, a point already 
noted by Fromang et al. (2002) in the context of proto- 
stellar discs. 

Eqns. - Q form a set of stiff coupled ordinary dif- 
ferential equations, and we use the Gear method for their 
solution at each point in the simulation domain and for 
each simulation time step. The major source of ionisation 
that we consider is X-rays from the central protostar, and 
our approach to calculating the ionisation rate is described 
below in Sect. |U 

2.2. Multifluid MHD equations 

As the mean-free path for collisions, and the ion gyro- 
radii, are very much smaller than the length scales we 
consider in our calculations, we adopt a multifluid MHD 



approach to incorporating chemical evolution of the gas 
during the dynamical evolution of our disc models. We 
use the shearing box representation of a local patch of 
the protostellar disc (Goldreich & Lynden-Bell 1965), in 
which the fluid is described using Cartesian coordinates 
(x, y, z). The origin of this coordinate system rotates 
with the local Keplerian angular velocity, fi, and the x 
coordinate represents the radial direction, the y coordi- 
nate the azimuthal direction, and z the vertical direction. 
The standard shearing box equations for MHD, including 
ohmic resistivity, are: 



dg 
di 

b\ 

~dt 
1. 



V-(gv) = 0. 



v • Vv 

-VP - — 
g gc 



-2flz x v + 3f2 2 a;x 
J x B - il 2 zi 



(5) 



dB 

~dt 



Vx (vxB — 77V x B) 



(6) 



The expressions above represent the continuity, momen- 
tum and induction equation, respectively, g represents the 
gas density, v the velocity, P the pressure, B the magnetic 
field, J the current density, and f2 is the local Keplerian 
angular velocity. The resistivity is denoted by 77. 

In our scheme the five species are treated as indivdual 
but tightly coupled fluids which move with the bulk veloc- 
ity v, and so we must solve a continuity equation for each 
of them. When combined with the possibility that the lo- 
cal abundance of species may change because of chemical 
evolution as well as advection, then the continuity equa- 
tion for each species i may be written: 



V- ( ft v) 



with i = l,. 



(7) 



at 

where m, is the particle mass for species i, Jj denotes the 
chemical reaction rate associated with the jth chemical 
reaction, while v^Js is the formation/destruction rate of 
the ith fluid component due to the jth chemical reaction. 
We assume an isothermal equation of state such that 



p = 4q 

and calculate the resistivity according to (Blaes 
1994) 



V : 



234 



T l/2_ 



(8) 
Balbus 

(9) 



The thermal structure of the disc model we use is de- 
scribed in the next section. 

2.3. Disc model 

The MHD simulations that we performed were calculated 
using a system of dimensionless variables, as is convenient 
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when performing shearing box simulations. In order to 
evolve the chemistry, however, whose reaction rates de- 
pend on local temperature and density, we need to ascribe 
physical units to these simulated quantities. We assume 
that the central star is of solar mass, and we adopt a disc 
model for which the surface density varies according to 
E = 1000 [R/l AU]~ 3/2 gem" 2 , and where the volume 
density varies with disc height as a Gaussian: 

e (z)ocexp|-i (|:) 2 }- (10) 

We assume that the disc aspect ratio has a constant value 
H/r = 0.05, and that the sound speed c s = HQ. The tem- 
perature is then related to the local sound speed according 
to c 2 = IZT/fi, where the mean molecular weight is as- 
sumed to be n — 2.33. The local physical parameters for 
each shearing box simulation are then fully specified once 
the orbital radius, magnetic field, chemical abundances 
and the local ionisation rate are defined. 

2.4. Numerical method 

We use the ZEUS finite difference MHD code (Stone & 
Norman 1992), modified to allow for the treatment of an 
arbitrary number of coupled fluids, to perform all our sim- 
ulations. Time step control was achieved using the stan- 
dard Courant condition, with the Courant number set to 
0.5. As the simulations use explicit time stepping, the time 
step size is also determined by the ohmic diffusion rate of 
the magnetic field. The chemical kinetic equations were 
evolved using the Gear method. 

2.5. Initial and boundary conditions 

The basic disc model used is described in Sect. 12.31 The ini- 
tial velocities of the gas in the simulations were taken to be 
the shearing box equilibrium values v = (0, —3/2 • ficc, 0), 
but with random fluctuations imposed with maximum 
amplitude equal to 10~ 3 of the sound speed. The initial 
magnetic field is a zero net flux vertical field given by 
B z = Bosm(2irx/H), where Bq is defined by the require- 
ment for the volume averaged plasma parameter (3 — 100. 

At the beginning of each shearing box simulation, 
the equilibrium particle concentrations 1^(7] of species 
Y G {Mg,Mg + ,H 2 ,H^} are taken as initial abundances. 
Note that we use different concentrations of magnesium 
for our models, and we simulate local patches of the disc 
at different radii from the central star. In each case we cal- 
culate the local equilibrium chemistry prior to initiating 
the MHD simulations. 

We use the same computational set-up as Fromang & 
Papaloizou (2006). The computational domain is given by 
[-H/2,H/2], [0,2%H], and [-3H,3H] in x, y, and z. We 
use a grid resolution of 32 x 100 x 192. Standard peri- 
odic boundary conditions apply in y and z, while periodic 
boundary conditions in shearing coordinates are used for 
x. For a detailed description of the shearing box set up 
and boundary conditions see Hawley, Gammie & Balbus 



(1995). Following Fromang & Papaloizou (2006) we intro- 
duced a vertical length scale H = 2.4 H which is used 
to prevent unphysical fluctuations due to the non vanish- 
ing vertical component of the gravitational force at the z 
boundary. By applying Eq. (30) of Fromang & Papaloizou 
(2006) we ensure that the vertical gravity acts on vertical 
length scales L < H only. 



3. Reaction-diffusion model 

In paper II we calculated the ionisation fraction for con- 
ventional a-disc models and examined the effect of turbu- 
lent mixing by modelling the diffusion of chemical species 
vertically through the disc. To recap: applying a one di- 
mensional reaction-diffusion model, we assumed that ver- 
tical mixing arises because of turbulent diffusion, and 
adopted the approximation T> = u t using the a prescrip- 
tion to calculate v t . Here v t is the (turbulent) kinematic 
viscosity that drives the radial diffusion of mass through 
the protostellar disc, and T> denotes the vertical diffusion 
coefficient. 

Instead of just mimicking the effects of turbulent mix- 
ing in this way, we now model the turbulent transport of 
chemical species by solving the corresponding non-ideal 
MHD equations in a three dimensional shearing box as 
discussed above. The mixing now is a direct outcome of 
the MHD turbulence. 

One purpose of this paper, however, is to examine 
whether or not the effects of turbulent mixing on the ioni- 
sation fraction described in paper II, can also be observed 
in shearing box simulations when MHD turbulence oper- 
ates. Hence, we calculated the ionistion fraction obtained 
for the kinetic model of Oppenheimer & Dalgarno by ap- 
plying the reaction-diffusion model at the corresponding 
cylindrical radius R, in order to aid a direct comparison 
between the results obtained for the reaction-diffusion and 
the shearing box model. We have good reason to make 
this comparison because Balbus & Papaloizou (1999) have 
shown that the mean flow dynamics of the MHD turbu- 
lence follows the a prescription. 

For the reaction-diffusion model we assume the same 
vertical density profile we use for the shearing box simu- 
lation at t — 0. The same applies for the gas temperature. 
We further adopt the approximation T> — v m with 

v m = a m c%/Q. (11) 

where a m is a dimensionlcss number used to quantify the 
efficiency of vertical mixing of chemical species by the tur- 
bulence. A range of a m values are used to examine how the 
reaction-diffusion model matches the simulations. Since 
we expect the largest gradients in the electron fraction to 
be in vertical (z) direction, we consider only vertical diffu- 
sion in the reaction-diffusion model. The rate of change of 
the molar density of the ith component of the fluid within 
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a given volume due to chemical reactions and diffusion 
caused by concentration gradients is 



dt 



d ( ^ d 
— I til) — x 
dz \ dz '' 



(12) 



where denotes the molar density of the ith component, 
T> the diffusion coefficient, x i — njn the fractional abun- 
dance of species i, and n = ^ rii. Jj denotes the chemical 
reaction rate associated with the jth chemical reaction, 
while v^Jj is the formation/destruction rate of the ith 
component due to the jth chemical reaction. A detailed 



discussion of the derivation of Eq. ( 12 ) is presented in pa- 
per II. 

The numerical method applied to solve the reaction- 
diffusion model is described in paper II, and the 
same boundary conditions apply. The boundaries are at 
[0, +3if ] such that the computational domain has size 
L z = 3H; the number of grid cells is L z /Az = 60, en- 
suring that the elements and charge are conserved to high 
accuracy. 

For a given metal elemental abundance, the reaction- 
diffusion models arc initiated with the equilibrium com- 
position obtained for T> — 0, exactly as they are for the 
shearing box simulations. 

4. X-ray ionisation rate 

As in paper I, we assumed that the ionisation of the disc 
material arises because of incident X-rays that originate in 
the corona of the central T Tauri star. We neglect contri- 
butions from Galactic cosmic rays as it remains uncertain 
whether they can penetrate into the inner disc regions we 
consider due to the stellar wind. The details for determin- 
ing the effective ionisation rate ( eS have been described 
in paper I. However, here we do not consider standard a- 
disc models as we did in that paper. Instead, we adopted 
a locally isothermal disc with the same Gaussian verti- 
cal profile for the mass density used for the local shearing 
box model at t = 0. When calculating the X-ray ionisation 
rate we assumed that this density profile did not vary with 
time, an assumption that is confirmed by the shearing box 
simulations which show that the density profile closely fol- 
lows its initial Gaussian profile throughout the nonlinear 
evolution of the MRI, (see top right panel of Fig. [3]) . 

We adopted values L x = 10 31 erg s _1 and k B T x = 
5 keV for the total X-ray luminosity L x and the plasma 
temperature T x , respectively. Compared with the val- 
ues applied in Ilgner & Nelson (2006a, b), (i.e. L x = 
10 30 erg s _1 and k B T x = 3 keV), the X-ray source consid- 
ered here is more both harder and more luminous in order 
to increase the ionisation fraction above a thcshold which 
makes the shearing box simulations feasible (a low ioni- 
sation fraction leads to a small time step size). However, 
the new values are still consistent with the observational 
constraints (e.g. Favata et al. 2005; Wolk et al. 2005). 

The effective ionisation rate (per hydrogen nucleus) 
C eff is approximated by Eq. (3) in paper I. Calculating the 




Fig. 1. The effective X-ray ionisation rate £ eff per hydro- 
gen nucleus. The contour lines refer to values of £ cff : 10~ 15 , 
_ i 

s 




z/H 



Fig. 2. Vertical profiles of the magnetic Reynolds numbers 
Re m . The profile drawn with the solid line refers to the 
profile used by Fleming & Stone (2003) in their model 
with a small dead zone. The two dashed lines correspond 
to profiles used by Turner et al. (2007) in their models 
Fl and F56. The profile drawn with the dotted line refers 
to our model 1, calculated at R = 10 AU and a metal 
abundance of x M = 5 x 10 -11 . Details are discussed in 
the text. 



X-ray optical depth r x along the line of sight between the 
X-ray source and the point in question, we derived the 
effective ionisation rate shown in Fig.[T] In particular, we 
applied the same data range in order to aid direct com- 
parison with the ionisation rates of our previous a-disc 
models (e.g. compare with Fig. 6 of paper I). The effec- 
tive ionisation rate shown in Fig. [l] is higher because of 
the brighter (by one order of magnitude) and harder (more 
penetrating) X-ray source applied. 

5. Previous simulations of dead zones 

As has been well documented in the literature, there re- 
main questions about the applicability of the MRI to pro- 
tostellar discs because of their high densities and low tern- 
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peratures, which lead to low levels of ionisation (e.g. Blaes 
& Balbus 1994; Gammie 1996). The first fully non lin- 
ear study of the MRI including the effects of resistivity 
were performed by Fleming et al. (2000), who employed 
shearing box simulations to examine the conditions under 
which fully developed turbulence could be sustained. They 
showed that the important quantity that determines the 
outcome is the magnetic Reynolds number, Re m , defined 
by: 

Rc m = ^ (13) 

n 

where the resistivity, T), is defined by Eq. (|9]). Simulations 
showed that turbulence could only be sustained if the 
Reynolds number was greater than a critical value, Re™ 1 *, 
where this critical value depends on the magnetic field 
topology. For zero net flux fields Re™* ~ 10 4 , whereas 
for net vertical fields Re™* ~ 100. The implication of this 
study is that protostellar discs, in which the magnetic field 
is an internally generated zero net flux field, will sustain 
turbulence in their surface regions where the ionisation de- 
gree causes Re m ^ 10 4 , but will remain in a near-laminar 
state in regions near the midplanc, as envisaged by the 
layered disc model of Gammie (1996). 

Shearing box simulations of stratified protostellar 
discs, with resistivity varying with height, were presented 
by Fleming & Stone (2003). Calculations were presented 
with different vertical resistivity and magnetic Reynolds 
number profiles, and it was shown that layered accretion 
resulted when the magnetic Reynolds number satisfied 
Re m > Re m * in the surface layers, with Re m < Re™* in 
the midplane regions. Their results also showed that a low 
Reynolds stress could be sustained in the dead zones due 
to the penetration of sound waves excited in the overly- 
ing active regions. The magnetic Reynolds number profile 
assumed in the simulation for a "small dead zone" from 
Fleming & Stone (2003) is shown by the solid line in Fig. [2] 
The magnetic Reynolds number varies from 1000 at the 
midplane, to 9.23 x 10 5 , 1.31 x 10 7 , and 1.62 x 10 s at 
z/H = 1, z/H = 2, and z/H = 3, respectively, resulting 
in a simulated disc with a dead zone whose vertical height 
is < H. 

In a recent paper, Turner et al. (2007) have presented 
a study of dead zones which included shearing box simu- 
lations of vertically stratified discs with resistivity varying 
as a function of height, and also a multifluid simulation in 
which the resistivity was able to change locally because of 
chemical evolution of the gas. In this study, Turner et al. 
(2007) employed the reaction network given by Eqns. (JIJ 
Q in order to calculate the ionisation fraction and re- 
sistivity. For those runs in which the resistivity was kept 
constant in time, the initial resistivity profile was obtained 
from the equilibrium solutions to Eqns. - Q- The run 
in which resistivity varied in time and space employed a 
multifluid approach, similar to that described in Sect. [2] 
Ionisation was assumed to be due to cosmic rays, and the 
underlying disc model was the minimum mass solar neb- 
ula model of Hayashi (1981). In Fig. [2] we present two of 



the resistivity profiles employed by Turner et al. (2007) 
corresponding to their runs Fl and F56. In run Fl the re- 
sistivity was a fixed function of height and corresponded to 
a radial location R = 1 AU with the gas-phase abundance 
of magnesium equal to the solar abundance (3.39 x 10~ 5 
magnesium atoms per hydrogen nucleus, corresponding to 
x[Mg] ~ 6.8 x 10~ 5 in our units). This run led to a lay- 
ered accretion flow with active surface layers and midplanc 
dead zone, as expected from the steep resistivity profile. 
It was shown that the boundary between dead and ac- 
tive zones is well described by the criterion that MHD 
turbulence is sustained by the MRI if the Lundquist num- 
ber Lu = v\/(£lrj) > 1, where va is the Alfven speed. 
Run F56 had a resistivity profile corresponding to a disc 
at 5 AU with gas phase magnesium abundance equal to 
10~ 6 below solar abundance (corresponding closely to our 
model with x[Mg] = 5 x 10~ n ). This model led to a fully 
turbulent disc, as expected from the flat resistivity profile 
obtained because the disc provides less shielding of cosmic 
rays at 5 AU. 

The multifluid model presented by Turner et al. (2007), 
in which the chemistry was evolved simultaneously with 
the dynamics, led to an interesting and somewhat unex- 
pected result. This disc model corresponded to the radial 
position R = 1 AU in the disc, and assumed a gas-phase 
magnesium abundance equal to the solar value. During 
the early phase of the model, the disc showed the expected 
layered structure with dead zone near the midplane, and 
active zone near the disc surface. After about 60 orbits 
the situation changed after a period of more intense mix- 
ing caused by enhanced turbulent activity. The recombi- 
nation time then exceeded the mixing time, allowing free 
electrons to mix toward the midplane, where net radial 
and azimuthal fields had built up due to field of the op- 
posite polarity advecting through the vertical boundaries. 
The presence of these net magnetic fields leads to an en- 
hance magnetic stress that partially enlivens the dead zone 
during periods when the ionisation fraction has been in- 
creased. 

In this paper, we assume that the ionisation of the 
disc material arises because of X-rays that originate in 
the corona of the central T Tauri star. We neglect contri- 
butions from Galactic cosmic rays, as it remains an open 
question whether or not they can to penetrate into disc 
regions we consider. The X-ray ionisation rate decreases 
with cylindrical radius R, because the X-ray optical depth 
along the line of sight increases as one moves out into the 
disc. 

Our use of strictly periodic boundary conditions in the 
vertical direction, along with an initial magnetic field that 
has zero net flux, means that the net flux remains zero 
throughout the simulations. The expectation then is that 
turbulence will be sustained only in those regions where 
the magnetic Reynolds number Re m <; 10 4 , and we do 
not expect that large scale net-flux magnetic fields will be 
able to accumulate in the dead zones of our simulations. 
As such we do not expect to observe the behaviour shown 
by simulation VI of Turner et al. (2007). We now present 
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the results from a series of systematic experiments which 
examine the effects of chemical evolution and turbulent 
mixing on the evolution of the MRI. 

6. Simulation results 

In this section we present the results of our multifiuid 
MHD simulations which examine the role of turbulent 
mixing on the structure of dead zones. The primary aim 
of these simulations is to demonstrate that there exists a 
region of parameter space in which turbulent mixing can 
enliven a dead zone which is otherwise predicted to exist 
in models that neglect turbulent transport. A further aim 
is to demonstrate good agreement between MHD simu- 
lations and the predictions of a simple reaction-diffusion 
scheme. This latter issue is addressed in Sect. 16.31 

When discussing the results of our simulations we 
will often refer to certain averaged quantities. We use 
the same averaging procedures presented in recent pub- 
lications studying vertically stratified disc models (e.g. 
Stone et al. 1996; Flemming & Stone 2003; Fromang k 
Papaloizou 2006). The azimuthal and radial average of 
quantity f(r, i) at a given time t is defined by 



f*om) = (/)* 



JJf(x',y',z,t) dx> dy' 



JJ dx' dy' 
and the volume average is given by 

JJJf(x',y',z',t) dx' dy' dz' 



F*'{t) = U) = 



JJJdx' dy' dz 



where the symbols (/)* and (/)**, denote the correspond- 
ing averaging procedures. The time-averaged values are 
denoted by (/)* and (/)**. 

A measure of the effective shear stress generated by the 
turbulence is given by the parameter a, which has contri- 
butions from both the Reynolds and Maxwell stresses: 



" - "Rey "T 

where 

a Rcy (x, z) 

aMax(x, z) 



'Max 



pRey , 

J£ §- = -wq( v x - (%)) K - (%)) 

M) 



T^Max 
± r<S> 

P 



-(B x B y ) 
4irP 



(14) 

(15) 
(16) 



and the azimuthal average (over y) is denoted by angled 
brackets, and P$ is the initial midplane pressure. We reg- 
ularly use the volume and time averaged values of a when 
discussing the simulation results below. In the following 
simulations we consider three different treatments of the 
resistivity and chemistry and we refer to these models as 
modell, model2 and model3. We describe each of these 
below, and each are summarised in table [2] 
1) modell : This model assumes a static resistivity profile 
which varies with height. The resistivity profile at t = 
is obtained using the equilibrium solution of the kinetic 
model presented in Eqns. Q - Q. During the simula- 
tions the local resistivity values are kept fixed. This model 



Table 2. List of models considered. Note that the col- 
umn "recombination process included" specifies whether 
or not recombination of free electrons occurs along with 
dynamical evolution. 



model 


resistivity 


recombination process included 


modell 


d t r](t, z) = 


no 


model2 


d t v(t, z) # 


yes 


mode!3 




no 



corresponds to the single fluid models of Fleming & Stone 
(2003), and the runs Fl, F52, F56, F58 of Turner et al. 
(2007). 

2) mode!2 : This is a multifiuid model in which the resistiv- 
ity varies in both time and space. The chemical reaction 
network given by Eqns. - Q is solved simultaneously 
with the dynamical evolution. The resistivity profile at 
t = is obtained from the equilibrium solution of the ki- 
netic model. 

3) mode!3 : This is a multifiuid model which has a resistiv- 
ity profile which varies in time and space. The resistivity 
profile at t = is obtained using the equilibrium solution 
of the kinetic model presented in Eqns. ([TJ - Q. For t > 0, 
however, the recombination of free electrons with ions is 
switched off, as is further ionisation of neutral species by 
X-rays. Local changes in the resistivity are due only to 
the turbulent mixing of ions. This model is equivalent to 
one in which the initial values of resistivity are conserved 
on fluid elements by being advected with the flow. 

We now present our simulation results in detail. We be- 
gin by highlighting simulations which show that turbulent 
mixing can remove the dead zone, before examining disc 
evolution at different radii and with different magnesium 
abundances. 



6.1. Dead-zone removal through turbulent mixing 

In this subsection we present a suite of models which 
demonstrate that turbulent mixing and continuing chemi- 
cal evolution of the gas are able to enliven a dead zone. All 
three models presented in this subsection correspond to a 
radial location in the disc R = 10 AU, and have a gas- 
phase magnesium abundance x Mg = 5 x 10 -11 . Because of 
the identical initial conditions used for these three models, 
we can identify the specific effects of the chemistry on the 
evolution of the MRI . 

modell 

The resistivity in this model is calculated from the equilib- 
rium electron abundance predicted by the chemical model 
presented in Sect. [2] and is held constant throughout the 
simulation. The variation of magnetic Reynolds number 
with height is shown in Fig. [2j which shows that the re- 
sistivity profile is intermediate between that applied by 
Fleming & Stone (2003) in their "small dead zone" model, 
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z/H 



Fig. 3. model l/model2 - Time averaged vertical profiles of the horizontally averaged density, normalised Maxwell and 
Reynolds stresses T R<i ,/P , and plasma parameter (3. The profiles on the left refer to the results obtained with the model 
assuming a fixed resistivity (modell). Results obtained with model2 are shown in the right panels. Both shearing box 



models are calculated for R = 10 AU and 



'Mg 



5 x 10 11 . a Max (solid line) and ct Rey (dashed - dotted line) refer to a 



time average taken over [100,200] orbits (for modell) and [0,100] orbits for model2. For all the other quantities, the time 
averages are taken over 10 orbit intervals, starting from t = (dashed line). For modell the solid line denotes averages 
takenbetween t = [140, 150] orbits, and for model2 the solid lines represented averages taken between t = [90, 100] 
orbits. The dotted lines refer to time averages taken over t — [0, 10], [10, 20], • • • , [80, 90] for model2, • • • , [130, 140] for 
modell. 



and the model F56 presented by Turner et al. (2007). 
Specifically Re m takes the following values: 1.4 x 10 3 at 
the disc midplane; 6.0 x 10 3 at z/H = 1; 1.4 x 10 4 at 
z/H = 2; 1.3 x 10 5 at z/H = 3. 



In agreeement with our expectations, this simulation 
resulted in a disc with well-defined dead and active zones, 
with the boundary between these it z/H w 1 

where Re m ~ 6 x 10 3 and Lu ~ 1. This dead zone is 
larger than that obtained by Fleming & Stone (2003) in 
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Fig. 4. model l/model2 - Time averaged vertical profiles of the horizontally averaged kinetic energy, 0.5gv 2 / Pq, and 
the magnetic energy B 2 /(8ttPq). The profiles on the left panel show results obtained assuming a fixed resistivity 
profile (modell), while the results obtained with model2 are shown in the right panels. Both shearing box models are 
calculated for R = 10 AU and x Mg = 5 x 1CP 11 . The time average is taken over 10 orbit intervals, starting from t = 
(dashed line) toward the interval t — [90, 100] (solid line, model2) and t = [140, 150] (solid line, modell), respectively. 
The dotted lines refer to time averages taken over t = [0, 10], [10, 20], • • • , [80, 90] ([130, 140]) orbits. 



their model whose resistivity profile is shown in Fig. [2] for 
which the dead zone was confined to z/H < 0.4. 

The time and volume averaged sum of the Maxwell 
and Reynolds stresses was found to be a** — 4.89 x 10 -3 
when the time average was taken over the interval [20,100] 
orbits. When taken over an interval [100,200] orbits, the 
stresses were found to decrease to a** = 8.49 x 10 -4 . This 
occurs because the stresses are higher during the develop- 
ment of the non linear stage of evolution early on in the 
simulation. We note that a fully active disc is expected to 
have a value a** 10 -2 (see later). 

The vertical profiles for the horizontally averaged den- 
sity q, the a values associated with the Maxwell and 
Reynolds stresses, the plasma parameter (3, the kinetic 
and magnetic energy are shown in the left hand panels 
of Figs. [3] and [4] The time averages are taken over 10 or- 
bit intervals, starting from t = (dashed line) towards 
[140, 150] (solid line). The dotted lines refer to profiles av- 
eraged over t = [0, 10], [10, 20], • • • [130, 140]. For each of 
the profiles, the same qualitative behaviour reported in 
Fleming & Stone (2003) is observed: 
(i) The vertical density profile remains unchanged 



throughout the nonlinear evolution of the MRI. 

(ii) A significant decline in the magnetic field energy to- 
wards the disc midplane is observed as compared to sur- 
face regions; j3 at z/H — is between 2 and 3 orders of 
magnitude greater than in the active zones. 

(iii) Dominance of the Reynolds stress over the Maxwell 
stress is observed in dead zones, due to the penetration of 
sound waves emitted in the overlying active zones, while 
the Maxwell stress is the dominant mechanism by which 
transport in active zones occurs. The transition occurs at 
z/H ~ 1 where the cross-over between active and dead 
zones occurs. 

mode!2 

In this model the initial free electron abundance and re- 
sistivity profile is calculated from the equilibrium solution 
to Eqns. Q - Q- The full set of multifluid equations, and 
the chemical network, are evolved together so that the lo- 
cal resistivity can change through turbulent transport of 
ions and chemical evolution (recombination/ionisation) of 
the gas. For this model we estimate that the turbulent 
mixing time corresponding to a ~ 0.01 is shorter than the 
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Fig. 5. model3: - Time averaged vertical profiles of the horizontally averaged value of normalised Maxwell and Reynolds 
stresses T R $/P , the plasma parameter (3, the kinetic energy 0.5gv 2 /P , an d the magnetic energy B 2 /(8nP ). As in 
the previous figures, the shearing box model is calculated for R = 10 AU and x Mg — 5 x 10~ n . a Max (solid line) and 
a Rcy (dashed - dotted line) refer to a time average taken over [0,100] orbits. For all the other quantities, the time 
average is taken over 10 orbit intervals, starting from t — (dashed line) towards [90,100] (solid line). The dotted 
lines refer to time averages taken over the intervals t — [0, 10], [10, 20], • • • , [80, 90] orbits. 



recombination time, such that chemical mixing should en- 
liven the dead zone. This is indeed what we find, as the 
simulation results in a turbulent flow which fills the entire 
volume of the disc, with a** = 1.31 x 10 -2 , where the time 
average was performed in the interval [20,100] orbits. We 
plot the vertical profiles of the various physical quantities 
that we have already described for modell in Figs. [3] and 
[4j Comparing the figures for modell and model2 we can 
make the following observations: 

(i) The mean vertical density profile remains approxi- 
mately constant in both models. 

(ii) Whereas the Reynolds stress dominates near the mid- 
plane in modell and the Maxwell stress dominates in the 
disc surface layers, we see that the Maxwell stress is dom- 
inant throughout in model2. 

(iii) The plasma (3 parameter is found to become very high 
(/3 ~ 10 5 ) in the midplane of the disc in modell as the 
magnetic field strength there becomes very low, whereas 
it remains in the range 10 < (3 < 10 3 for \z/H\ < 2 for 
model2, indicative of a fully active disc. 

(iv) The kinetic energy throughout the disc, but espe- 
cially near the midplane, is much higher in model2 than 



in modell as the turbulent velocity field is driven by the 
MRI. 

(v) The magnetic energy near the disc midplane for 
mode 12 is more than two orders of magnitude greater than 
for modell due to the continuing dynamo action associ- 
ated with the MRI. 

We conclude that model2 shows unambiguously that the 
dead zone in the disc can be removed by turbulent mixing 
and continuing chemical evolution under circumstances 
where the local mixing time scale is smaller than the re- 
combination time, in basic agreement with the prediction 
of the reaction-diffusion model presented in paper II. 

mode!3 

At t — the initial resistivity profile is set up in the same 
way as described for modell and model2. For t > 0, the 
local resistivity is updated after every MHD time step be- 
cause of the transport of free electrons and ions only. Due 
to the inhibition of recombination and ionisation in this 
model, free electrons diffuse and cause the resistivity to 
become homogeneous on long time scales. In the presence 
of sufficient numbers of free electrons in the initial ionisa- 
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tion state of the disc, we expect that mixing will lead to a 
fully active disc, and indeed this is what we find. Instead 
of the two-layer structure obtained using model 1 above, 
with dead and active zones, the MHD turbulence now fills 
the full vertical extent of the disc. For t > 40 orbits, we 
observe a quasi steady state characterised by small fluctu- 
ations around the mean value a** = 5.75xlCP 3 , where the 
time average was taken in the interval [20,100] orbits. The 
vertical profiles of various quantities are shown in Fig. [5] 
Compared with the results obtained for model 1, we see 
that the volume averaged turbulent stresses in each case 
are very similar, even though modell had a dead zone. 
The reason for this is that the resistivity in model3 is 
higher near the disc surface because of mixing, and so 
reduces the strength of the turbulence there. The subse- 
quent enlivening of the midplane in model3 does not lead 
to a substantial increase in the overall stress because the 
now-uniform resistivity is sufficient to damp the strength 
of the turbulence compared to its state in an ideal MHD 
calculation. 

Comparing model3 with model2 we see that higher 
stresses and more vigorous turbulence are generated by 
model2. This is because model3 generates a disc without 
a dead zone, but in which the resisivity is higher than for 
model2, such that the strength of the resulting turbulence 
is suppressed somewhat. The results for model3 show that 
mixing the initial free electron population throughout the 
disc can cause the dead zone to disappear, but that allow- 
ing the chemical evolution of the disc to continue during 
the turbulent mixing leads to a more active disc. This is 
because the continuing ionisation of species near the disc 
surface, followed by mixing toward the midplane, produces 
a higher ionisation fraction overall. 

6.2. Results as function of xyig and orbital radius 

A primary motivation for this paper was the indication 
in paper II that dead zones could be enlivened by a com- 
bination of turbulent mixing and sufficient abundance of 
gas-phase magnesium atoms. In that paper we presented 
calculations of the ionisation fraction in standard a-discs 
using reaction-diffusion models. The main results were 
that turbulent mixing could only change the structure of 
a dead zone if: (i) the abundance of magnesium was suf- 
ficient (so as to increase the recombination time); (ii) one 
was considering locations further out in the disc where the 
lower temperatures and densities increase the recombina- 
tion time relative to the local turbulent transport time. 

The purpose of the simulations presented in the fol- 
lowing subsections are to examine how turbulent mixing 
affects dead zone structure as a function of magnesium 
abundance and radial position in the disc, as a test of 
the predictions contained in paper II. We also perform a 
detailed comparison between some of our MHD simula- 
tions and the predictions of the reaction-diffusion model. 
These simulations solve the full set of multifluid equa- 
tions in combination with the chemical model described in 



Sect. [2] We first present shearing box simulations of discs 
at various locations between 1 and 10 AU with gas-phase 
magnesium abundance equal to zero, before considering 
a similar set of models with gas-phase magnesium abun- 
dance xm s = 5 x 10~ n , which is about 10 -6 of the solar 
value. 

6.2.1. Models with x Mg = 

We begin our discussion by first examining the dead 
zone structure after saturation obtained for model2 with 
XMg — 0. We examine the disc evolution at radii R = 1, 
3, 5, 7 and 10 AU. In basic agreement with the results 
obtained in paper II, the dead zone structures obtained 
when magnesium is absent are very similar for all radii 
considered. In Fig. [6] we present space-time plots of the 
horizontally averaged value of a, and it is clear that the 
disc sustains a two-zone structure for all radii and all time, 
which consists of a large dead zone which extends from 
the midplane to \z/H\ ~ 2 where the magnetic Reynolds 
number Re m > 4000 and the Lundquist number Lu <; 1 
(shown by the solid lines in Fig. [fj|. Across the region 
bounded by 2 < \z/H\ < 2.3 we find that the horizontally 
averaged a varies by more than two orders of magnitude, 
and this is maintained for the duration of the simulation 
(100 orbits). 

6.2.2. Models with x Mg = 5 x 1CT 11 

We now consider the dead zone structures obtained for a 
magnesium abundance of Xy[g = 

5 x lO" 11 at radii R = 
1, 3, 5, 7, and 10 AU. Space-time plots of the hoizontally 
averaged value of a are shown in Fig. [7} At R = 1 AU and 
3 AU, the location of the boundary separating the dead 
from the active layer matches very well the structure ob- 
tained for ir Mg = 0, since the height at which Re m begins 
to exceed 4000 is \z/H\ ~ 2 (which is also the height at 
which the Lundquist number Lu > 1, shown by the solid 
lines in Fig. [7]). Significant changes in the dead zone struc- 
ture start to become evident at R > 5 AU. At this radius 
the recombination time is similar to the turbulent mixing 
time, allowing the resistivity to be reduced there through 
the transport of free electrons into the dead zone. Over 
longer time scales we see that the dead zone that is estab- 
lished early on in the simulation starts to diminish, and 
between 80-100 orbits there is evidence that the dead zone 
size has decreased significantly. Nonetheless, at 100 orbits 
we find that there remains a region in the vicinity of the 
midplane that retains an average values of a ~ 10~ 4 . 

By contrast, the space-time plots of the horizontally 
averaged value of a at R = 7 and 10 AU show that the 
two-layer structure has completely disappeared, because 
mixing reduces the resistivity, and hence increases the 
magnetic Reynolds number to values <; 4000. Turbulence 
now fills the entire computational domain which is con- 
firmed by the time and volume averaged values of a**: 
1.07 x 10" 2 for R = 7 AU and 1.31 x 10~ 2 for R = 10 AU. 
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Our results indicate that dead zones can be reduced 
or removed altogether through turbulent mixing. The 
criteria for sustaining MHD turbulence have already been 
discussed in the context of reaction-diffusion models 
in paper II, and also apply to the shearing box models 
considered here. These criteria are: 

(i) There are sufficient metal atoms available in the gas 
phase so that recombination between metal ions and 
electrons becomes the dominant process by which the 
local ionisation fraction is determined. 

(ii) The turbulent mixing time scale is shorter than the 
dominant recombination time on which free electrons are 
removed. 

We now compare in detail the results of some of our 
MHD simulations with those obtained using the reaction- 
diffusion model presented in paper II and discussed in 
Sect. ED 

6.3. Comparing shearing box simulations and 
reaction-diffusion models 

Encouraged by the good qualitative agreement obtained 
between the MHD simulations and the results presented 
in paper II, we now examine in detail the level of agree- 
ment between the simulations and the reaction-diffusion 
model. As mentioned in Sect. [3] we assume that the diffu- 
sion coefficient, T>, which governs the rate at which chem- 
ical species mix vertically in the disc, is equal to an ef- 
fective kinematic viscosity generated by the turbulence 
v m = oi m c 2 s /Q,, where ot m refers to a dimensionless pa- 
rameter that measures the rate of vertical mixing (not to 
be confused with the value of a associated with the ra- 
dial transport of angular momentum). When solving the 
reaction-diffusion equations we use a range of ot m values 
to obtain different solutions, which we then compare with 
the results of the MHD simulations. Previous work on the 
vertical mixing of dust grains and molecules (Carballido 
et al. 2005; Johansen & Klahr 2005; Turner et al. 2006; 
Fromang & Papaloizou 2006) suggests that the ratio of the 
rate of angular momentum transport to the rate of trans- 
port of chemical species by the turbulence should lie in the 
range v/v m — 1 _ 3, and we examine how well our best-fit 
value of v rn agrees with this expectation. Furthermore, the 
results of Turner et al. (2006) and Fromang & Papaloizou 
(2006) show that the diffusion coefficient is not constant 
with height above the midplane, but increases with height 
because the turbulent velocities increase in proportion to 
the Alfven speed. Although we consider only a constant 
value of T> for each reaction-diffusion model that we run, 
we examine the quality of the best -fit that we obtain, and 
quantify this by stating the error obtained in predicting 
the magnetic Reynolds number (which is a proxy for the 
free electron abundance) at the disc midplane and surface. 
In general we find that using a uniform diffusion coefficient 
leads to a slight overestimate of the mixing rate near the 
midplane, and a slight underestimate near the surface. 



Overall we find good agreement between the reaction- 
diffusion model and the MHD simulation when the value 
of ce m adopted in the former is between a factor of 1-2 
times lowerer than the time and volume averaged value of 
a obtained in the latter. The a values corresponding to 
each MHD simulation are listed in table [3] The best-fit 
values of a m are listed in table [3] along with the Schmidt 
number which measures a/a m . 

We have compared the results of models located at 
disc radii R = 5, 7 and 10 AU, and with magnesium 
abundances xug, = and XMg = 5 x 10~ n . In each 
case we examine how the vertical profile of the magnetic 
Reynolds number Re m evolves with time, and in the case 
of the reaction-diffusion model we calculate the equilib- 
rium value of Re m from the underlying ionisation fraction. 
In each case we initiate the calculation assuming that the 
initial chemical abundance profile is equal to the equilib- 
rium state in the absence of mixing. For each disc radius, 
we plot the profile of the evolving Re m below. In the left 
panels we present the results from the MHD simulations, 
and in the right panels the equilibrium profile from the 
corresponding reaction-diffusion model assuming different 
values of a m . The values of Re m plotted for shearing box 
simulations are horizontal and time averages, where each 
time average was performed over 10 orbit intervals. We 
plot the initial value of Re m using the dashed line, and 
subsequent values are plotted using dotted lines starting 
at t = [0,10] and moving up to t = [80,90]. The final 
values at t = [90, 100] are shown using the solid line. 



Table 3. Time and volume averaged values of a** ob- 
tained for the models described in the paper. The first 
column gives the model label and the cylindrical radius 
considered. The values of a** are listed in the 2 nd and 
3 rd column assuming a elemental metal abundance XMg of 
and 5 x 10 -11 , respectively. Apart from the value de- 
noted with the symbol f , the time average was taken over 
[20,100] orbits, while ()^ refers to the time average taken 
over [100,200] orbits. 



a**(x Mg = 0) a"(% g = 5 • 10 xl ) 



modell 



at 10 AU 






4.89 


10~ 3 








(8.49 • 


w-y 


model3 










at 10 AU 






5.75 


10~ 3 


model2 










at 1 AU 


2.67 


10~ 4 


2.51 


10" 4 


at 3 AU 


3.01 


10" 4 


3.46 


10" 4 


at 5 AU 


2.84 


10" 4 


1.11 


10" 3 


at 7 AU 


5.44 


10~ 4 


1.07 


10" 2 


at 10 AU 


5.10 
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1.31 


10" 2 
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Table 4. Values of the diffusion coefficients applied to 
the reaction-diffusion model which best matches the cor- 
responding MHD results for model2 with XMg = 5 • 10~ n . 
The Schmidt number S c is listed in the 3rd column refer- 
ing to time and volume averaged values of a** between 
[20,100] orbits. 



R [AU] 


a 


111 


Sc 


5 


5.89- 


1(T 4 


1.88 


7 


8.67- 


1(T 3 


1.23 


10 


9.28- 


icr 3 


1.41 



6.3.1. Results at 5 AU 

We first discuss the results for R = 5 AU which are shown 
in Fig. [8] Note that in the upper left panel of Fig. [8] 
the magnetic Reynolds number shows a well defined min- 
imum, and this arises because this MHD simulation was 
performed with a ceiling being adopted for the resistivity 
in the induction equation, corresponding to a minimum 
value of Re m = 20. This was done to ensure that the time 
step constraint arising from the diffusive term in the in- 
duction equation was not too severe. The minimum value 
of Re m that we calculated from the electron fraction dur- 
ing this simulation was Re m = 10.007. 

The upper left panel of Fig. [8] simply shows that the 
model at 5 AU with XMg = maintains a significant dead 
zone with \z/H\ < 2 (where Re m ~ 4000) throughout, 
and the magnetic Reynolds number does not change from 
its initial value. The upper right hand panel shows that 
the reaction-diffusion equation agrees with this, as the 
single line plotted is actually three lines overplotted cor- 
responding to a m — 0, 10 -3 and 10~ 2 . In the absence of 
magnesium, the recombination rate is simply too high to 
allow mixing to modify the dead zone for a m values in 
this range. 

The lower left panel shows the evolution of Re m with 
XM g = 5x 10~ n . It is clear that the Re m profile in this case 
is non stationary near the midplane, even after 100 orbits, 
and this appears to be because this particular model is 
one which maintains a dead zone throughout the run, but 
whose parameters are close to those which would allow 
mixing to remove the dead zone. Episodic increases and 
decreases in turbulent activity modify the ionisation state 
near the midplane, causing the Re m values to oscillate 
about a value close to 1000. 

The lower right panel shows the results from the 
reaction-diffusion model agree quite well with the mean 
Rc m profile from the MHD run, in particular when a m = 
10~ 3 . Inspection of table [3] shows that mean value of a 
obtained from the MHD run was 1.52 x 10~ 3 . Table [I] 
shows that the best fit reaction-diffusion model has a 
value of a m — 5.89 x 10~ 4 (such that the Schmidt number 
equals 1.88). The error in the predicted value of Re m at 
the midplane was 3 %, while the error at the disc surface 
was ~ 15%, showing good overall agreement even when 



a uniform diffusion coefficient is adopted in the reaction 
diffusion models. 

6.3.2. Results at 7 AU 

The vertical profiles of Re m obtained at R — 7 AU are 
shown in Fig. [9] The upper left and right panels show re- 
sults from the MHD simulation and reaction-diffusion cal- 
culation, respectively for magnesium abundance XMg = 0. 
Once again we see that mixing has no effect on the mag- 
netic Reynolds number profile in the absence of magne- 
sium, and a substantial dead zone is maintained through- 
out both calculations, which show excellent agreement. 

The lower left and right panels show models for which 
XMg = 5 x 10~ n . Here we see that there is very significant 
change in the magnetic Reynolds number profile as tur- 
bulent mixing ensues. In the MHD simulation we see that 
the minimum value of Re m changes from 1000 to ~ 4000, 
which is high enough for the disc midplane to become ac- 
tive. The lower right panel shows good agreement with the 
MHD simulation when a — 10 -2 . We see from table[3]that 
the average value of a from the MHD run is 1.07 x 10~ 2 . 
Table [4] shows that the best fit reaction-diffusion model 
has a value of a m = 8.67 x 10 -3 (such that the Schmidt 
number equals 1.23). The error in the predicted value of 
Re m at the midplane was 2 %, while the error at the disc 
surface was ~ 15%. 

6.3.3. Results at 10 AU 

The profiles obtained at R = 10 AU are shown in Fig. [10] 
The upper panels are again in good agreement when 
XMg = 0, showing that mixing has no effect on the dead 
zone structure. The lower panels show that the Re m profile 
is changed significantly by mixing when 10" 11 , 
such that the dead zone is enlivened completely In the 
MHD simulation the minimum value of Re m changes from 
~ 1000 to ~ 6000, allowing the dead zone to become MRI- 
active and the disc to be turbulent throughout its height. 
Good agreement in the Re m profile is obtained using the 
reaction-diffusion model when a m ~ 10~ 2 , which as ex- 
pected is slightly lower than the value a = 1.31 x 10 -2 
listed in table [3] as arising from the MHD simulation. 
Table [4] shows that the actual best fit reaction-diffusion 
model has a value of a m = 9.28 x 10 -3 (such that the 
Schmidt number equals 1.41). The error in the predicted 
value of Re m at the midplane was < 1 %, while the er- 
ror at the disc surface was ~ 15%, which again illustrates 
the fact that reasonable agreement can be obtained when 
using a uniform diffusion coefficient. 

6.3.4. Reaction-diffusion results for model3 

We finally present a comparison between the MHD simu- 
lation for model3 and a corresponding reaction-diffusion 
model. To recap: model3 allows the free electrons and ions 
to diffuse, but does not include recombination or on-going 
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ionisation. The MHD simulation and reaction-diffusion 
model were initiated with the equilibrium chemical abun- 
dance for the case iM g = 5x 10~ n at R = 10 AU. The 
expectation is that turbulent mixing will cause the Re m 
profiles to change from their initial values to become uni- 



form. Inspection of Fig. 1 1 confirms that our models agree 
with this expectation. 

7. Conclusions 

We have presented the results from a scries of shearing 
box multifluid MHD simulations aimed at examining the 
evolution and structure of dead zones in protoplanetary 
discs, in the presence of turbulent transport of ions, on- 
going chemical evolution of the gas, and ionisation due to 
X-rays emitted by the central star. We have adopted a 
number of simplifying assumptions, including the absence 
of small grains whose presence in even modest numbers 
would lead to rapid removal of free electrons (Sano et al. 
2000; Ilgner & Nelson 2006a). As such, our results are 
likely to be most applicable to protostellar discs at a fairly 
late stage of evolution after grains have accumulated to 
form larger bodies. 

A primary objective of this work was to use MHD 
simulations to re examine the results of Ilgner & Nelson 
(2006b), who used a simple reaction-diffusion model to 
calculate the effects of turbulent mixing on dead zone 
structure. The simple model predicted that turbulent 
transport can be effective at enlivening a dead zone pro- 
vided that: (i) the abundance of gas-phase magnesium 
is sufficient; (ii) one considers regions of the disc beyond 
radii 5 AU where turbulent mixing occurs on a shorter 
time scale than recombination. The main conclusions of 
this paper are that full multifluid MHD simulations are 
in good agreement with these predictions. Models simu- 
lated at radii between 1-10 AU, and with no magnesium 
in the gas-phase, showed a two-layer structure consisting 
of an actively accreting zone near the disc surface, and a 
magnetically inactive region near the midplane. The ad- 
dition of gas-phase magnesium with fractional abundance 
XM g = 5 x 10~ n led to significant dead zones persisting 
for radii R < 5 AU, but models at 7 and 10 AU resulted 
in fully active discs without dead zones. The implications 
for protoplanetary discs is that at late times, when most 
of the small submicron sized dust grains have grown to 
much larger sizes, the dead zone beyond 5 AU may be en- 
livened because of turbulent transport of ions toward the 
midplane. Regions interior to 5 AU will, however, retain 
their dead zones. 

A further conclusion of our work is that detailed com- 
parison between the simple reaction-diffusion model and 
the MHD simulations leads to very good agreement in the 
vertical profiles of resistivity and magnetic Reynolds num- 
ber when an appropriate diffusion coefficient is chosen. 
Typically we find that the best-fit vertical diffusion coef- 
ficient corresponds to a ratio in the range 1-2 between the 
rate at which angular momentum is transported radially 
and the rate at which diffusion of chemical species occurs 



vertically. This result is consistent with those presented 
by Carballido et al. (2005), Johansen & Klahr (2005), 
Turner et al. (2006) and Fromang & Papaloizou (2006) 
who showed that turbulent diffusion of dust (and chemi- 
cal tracers) occurs on a slightly slower time scale than the 
transport of angular momentum, since it is driven through 
correlations in the perturbed flow velocities only. 

It has traditionally been assumed that MRI-driven 
MHD turbulence in discs can be sustained against the 
damping effects provided by ohmic resistivity if magnetic 
field diffusion over the characteristic wavelength of the in- 
stability occurs on a time scale longer than the growth 
time. Indeed Turner et al. (2007) showed that such a con- 
dition provides a good indicator of where the transition 
between active and dead zones in a disc will occur. They 
showed that the transition zone occurs where the Lunquist 
number Lu = v\j (rfit) ~ 1, and our simulation results are 
in good agreement with this. Recent work by Fromang & 
Papaloizou (2007) and Fromang et al. (2007), however, has 
shown that in the case of non stratified shearing box sim- 
ulations, turbulence is only sustained in discs where the 
magnetic Prandtl number P m = v jr\ > 1 (where v is the 
physical (molecular) viscosity), even when Lu > 1 in the 
initial state. The interpretation is that MRI-driven tur- 
bulence cascades the magnetic field down to the smallest 
scales available, by virtue of the turbulent velocity field 
twisting the field up. If the characteristic scale on which 
velocity fluctuations are damped is smaller than the resis- 
tive scale, then a zero net flux field will be dissipated and 
turbulence will die. Interestingly, the magnetic Prandtl 
number in protostellar discs is expected to be P m < 1, 
since resistivity is high and viscosity is low. Fromang & 
Papaloizou (2007) also show that the intrinsic numerical 
magnetic Prandtl number of the ZEUS code is > 1, at 
least for simulations with resolutions feasible on current 
computers. This suggests that the results in this paper, 
and those in other papers that have looked at dead zones, 
are modelling discs which only fullfil one of the necessary 
criteria for MHD turbulence to be sustained in a physical 
way, namely that Lu > 1. The condition for P m > 1 is 
satisfied because of the nature of numerical dissipation in 
the code. We note that the effects observed by Fromang 
& Papaloizou (2007) and Fromang et al. (2007) occur for 
the particular case of non stratified shearing box simula- 
tions, and that a mechanism for maintaining active MRI 
turbulence may be the generation of large scale magnetic 
field through magnetic buoyancy effects and field stretch- 
ing in vertically stratified discs, such as those we consider 
in this paper. Nonetheless, it is clearly necessary to ex- 
amine these issues by including the appropriate viscous as 
well as resistive transport coefficients in simulations, and 
we will do this in a future publication. 

There are two additional issues that we have not ad- 
dressed in this paper. The first is that scattering of X- 
rays toward the disc midplane may increase the ionisation 
rate in the disc by up to an order of magnitude (Igea & 
Glassgold 1999), and this can have an obvious effect on 
the structure of the dead zone. Although we have not un- 
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dertaken an extensive analysis of the effect of this, we have 
run a model at 1 AU with XMg = 5xlCP n with the X-ray 
luminosity increased by two orders of magnitude. We find 
only a small change in the dead zone structure in this case. 
We would expect in general that increases in the X-ray lu- 
minosity due to scattering will move the radial boundary 
of the dead zone inward slightly, but will not completely 
remove the dead zone. A final issue that we have not ad- 
dressed in this paper is that of X-ray flares. Observations 
of T Tauri stars by CHANDRA have shown that they emit 
regular outbursts of X-rays which may increase the X-ray 
luminosity by a few orders of magnitude, and also harden 
the X-ray spectrum (Favata et al 2005; Wolk et al. 2005). 
This issue was examined by Ilgner & Nelson (2006c) , who 
showed that the flaring could significantly modify dead 
zones in protoplanetary discs. We will address this issue 
in a future paper using multifluid MHD simulations with 
chemistry, so that both the effects of X-ray flaring and 
chemical mixing on dead zone structure can be examined. 
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Fig. 6. model2: - Space-time plots of the horizontally averaged value of a calculated at different radial positions R and 
a metal elemental abundance of a; M = 0. The cylindrical radius R varies from 1 AU (top panel) to 10 AU (bottom 
panel) passing through R — 3, 5, and 7 AU. The solid line drawn at the boundary between the live and dead zones 
indicates the position where the horizontally averaged Lundquist number Lu = 1 . 
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Fig. 7. model2: - Space-time plots of the horizontally averaged value of a calculated at different radial positions R 
and a metal elemental abundance of x Mg — 5 x 10 -11 . The cylindrical radius R varies from 1 AU (top panel) to 10 
AU (bottom panel) passing through R — 3, 5, and 7 AU. The solid line drawn at the boundary between the live and 
dead zones indicates the position where the horizontally averaged Lundquist number Lu = 1. 
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Fig. 8. model2: - Comparison between vertical profiles of the magnetic Reynolds number obtained for R = 5 AU 
and x Mg = and x M = 5 x 1CP 11 , respectively. Left panels show shearing box simulation results, and right panels 
reaction-diffusion model results. In the left panels, the time averaged vertical profiles of the magnetic Reynolds number 
Re* are shown by taking the time average over 10 orbit intervals, starting from t = (dashed line) towards [90, 100] 
(solid line). The other profiles (dotted lines) refer to the intervals t = [0, 10], [10, 20], • • • [80, 90] orbits. The resistivity 
profiles of the magnetic Reynolds number Re m shown in the right panels refer to the equilibrium profiles obtained 
with the reaction diffusion model assuming different values of a m . In particular for iMg — 0, the value of the magnetic 
Reynolds number is not affected by the actual value of a m , leading to identical profiles for a m = 0, a m = 10~ 3 , and 
a m = 10~ 2 . Note that both shearing box simulations and the reaction-diffusion models are initiated with the steady 
state profile obtained for a m = 0. 
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Fig. 9. model2: - Comparison between vertical profiles of the magnetic Reynolds number obtained for R = 7 AU 
and x Mg = and x M = 5 x 1CP 11 , respectively. Left panels show shearing box simulation results, and right panels 
reaction-diffusion model results. In the left panels, the time averaged vertical profiles of the magnetic Reynolds number 
Re* are shown by taking the time average over 10 orbit intervals, starting from t = (dashed line) towards [90, 100] 
(solid line). The other profiles (dotted lines) refer to the intervals t = [0, 10], [10, 20], • • • [80, 90] orbits. The resistivity 
profiles of the magnetic Reynolds number Re m shown in the right panels refer to the equilibrium profiles obtained 
with the reaction diffusion model assuming different values of a m . In particular for iMg — 0, the value of the magnetic 
Reynolds number is not affected by the actual value of a m , leading to identical profiles for a m = 0, a m = 10~ 3 , and 
a m = 10~ 2 . Note that both shearing box simulations and the reaction-diffusion models are initiated with the steady 
state profile obtained for a m = 0. 
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Fig. 10. model2: 



Comparison between vertical profiles of the magnetic Reynolds number obtained for R 
-li 



10 AU 



and x Mg = and x Mg = 5x 10 , respectively. Left panels show shearing box simulation results, and right panels 
reaction-diffusion model results. In the left panels, the time averaged vertical profiles of the magnetic Reynolds number 
Re* are shown by taking the time average over 10 orbit intervals, starting from t = (dashed line) towards [90, 100] 
(solid line). The other profiles (dotted lines) refer to the intervals t = [0, 10], [10, 20], • • • [80, 90] orbits. The resistivity 
profiles of the magnetic Reynolds number Re m shown in the right panels refer to the equilibrium profiles obtained 
with the reaction diffusion model assuming different values of a m . In particular for iMg — 0, the value of the magnetic 
Reynolds number is not affected by the actual value of a m , leading to identical profiles for a m = 0, a m = 10~ 3 , and 
a m = 10~ 2 . Note that both shearing box simulations and the reaction-diffusion models are initiated with the steady 
state profile obtained for a m = 0. 
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Fig. 11. model3: - Comparison between vertical profiles of the magnetic Reynolds number obtained for R = 10 AU 
and x Mg = 5 x 1CP 11 obtained by the shearing box simulation (left panel) and the reaction-diffusion model (right 
panel). The results refer to the simulation where recombination and ionisation are switched off (model3). In the left 
panel, the time averaged vertical profiles of the magnetic Reynolds number ReJ^ are shown by taking time averages 
over 10 orbit intervals, starting from t = (dashed line) towards [90, 100] (solid line). The other profiles (dotted lines) 
refer to t = [0, 10], [10,20], • • • [80,90]. The resistivity profiles of the magnetic Reynolds number Re m shown in the 
right panel refer to the equilibrium profiles obtained with the reaction-diffusion model assuming different values of 
a m . In fact, the equilibrium profile obtained for model3 is not affected by the actual value of a m > 0. Note that 
both shearing box simulation and the reaction-diffusion model are initiated with the equilibrium profile obtained for 



